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ABSTRACT 

In this paper we compare fully compressible (Meakin & Arnett 2006a, b) and anelastic (Kuhlen, 
Woosley, & Glatzmaier 2003) simulations of stellar oxygen shell burning. It is found that the two 
models are in agreement in terms of the velocity scale {vc ~ 10'' cm/s) and thermodynamic fluctuation 
amplitudes (e.g., p' / {p) ^ 2 x 10~^) in the convective flow. Large fluctuations ('^11%) arise in the 
compressible model, localized to the convective boundaries, and are due to internal waves excited in 
stable layers. Fluctuations on the several percent level are also present in the compressible model 
due to composition inhomogeneities from ongoing entrainment events at the convective boundaries. 
Comparable fluctuations (with amplitudes greater than ^1%) are absent in the anelastic simulation 
because they are due to physics not included in that model. We derive an analytic estimate for 
the expected density fluctuation amplitudes at convective boundaries by assuming that the pressure 
fluctuations due to internal waves at the boundary, p'^, balance the ram pressure of the convective 
motions, pv^. The predicted amplitudes agree well with the simulation data. The good agreement 
between the anelastic and the compressible solution within the convection zone and the agreement 
between the stable layer dynamics and analytic solutions to the non-radial wave equation indicate 
that the compressible hydrodynamic techniques used are robust for the simulated stellar convection 
model, even at the low Mach numbers found M ^ 0.01. 
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1. INTRODUCTION 

Oxygen burning (by ^^O -|-^^ O fusion) occurs in the 
precollapse stages of the evolution of massive stars. Neu- 
trino cooling speeds these stages to the extent that 
the evolutionary times scales are close enough to the 
sound travel time so that direct compressible numeri- 
cal hydrodynamics can be applied (Arnett 1994). The 
flrst detailed studies of this stage (Bazan & Arnett 
1998) were done in two-dimensional symmetry (2D) with 
PROMETHEUS (FryxeU, MiiUer, & Arnett 1989), a 
multi-fluid multidimensional compressible hydrodynam- 
ics code based on the Piecewise Parabolic Method (PPM) 
of Colella & Woodward (1984). They showed vig- 
orous convection, with significant density fluctuations 
(up to 8%) at the convective-nonconvective boundaries. 
These results were confirmed in detail in 2D with the 
VULCAN Arbitrary-Lagrangian-Eulerian (ALE) hydro- 
dynamics code (Livne 1993) by Asida & Arnett (2000). 
VULCAN is an entirely independent compressible hy- 
drodynamics code, so that these two sets of simulations 
only shared the initial model, the sonic time step lim- 
itation, and the 2D geometry. A new version of the 
PROMETHEUS code, PROMPI (which uses the Mes- 
sage Passing Interface for parallelism) , has extended the 
study to 3D. In all these compressible models except the 
earliest (Arnett 1994) the computational domain has in- 
cluded both the convective oxygen burning shell as well 
as two bounding stably stratified layers. 

Kuhlen, Woosley, & Glatzmaier (2003) investigated 
shell oxygen burning in 3D using an anelastic hydrody- 
namics code which filters out sound waves and linearizes 
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thermodynamic fluctuations around a background refer- 
ence state (Glatzmaier 1984; Cough 1969). In contrast to 
the fully compressible results above, Kuhlen, Woosley, & 
Glatzmaier (2003) found only small density and pressure 
contrasts, and subsonic flows which were well within the 
anelastic approximation (all thermodynamic contrasts 
less than 1 percent). The boundary conditions used were 
impermeable and stress-free and were placed within the 
convection zone so that convective overshoot could not 
be studied. In particular, the dynamic consequences of 
the neighboring non-convective shells, and their elastic 
response to convective fluctuations, were ignored. The 
formulation was single fluid, so that effects depending 
upon composition, i.e., mixtures of fuel and ashes, were 
not modeled. 

The applicability of both fully compressible and anelas- 
tic hydrodynamic methods have recently been challenged 
by developers of low-Mach number solvers (Almgren et 
al. 2006). The reliability of compressible codes has been 
questioned for low velocity flows due to possible viola- 
tions of elliptic constraints that arise in the evolution 
equations in the very small Mach number limit (e.g., 
Schneider et al. 1999). The limits for which compressible 
solvers remain robust in the astrophysical context, how- 
ever, has not been rigorously studied. Anelastic methods, 
on the other hand, enforce a divergence constraint on the 
velocity field which filters out sound waves, but are for- 
mulated assuming that thermodynamic fluctuations are 
small and only linear deviations from a background ref- 
erence state are retained. Therefore, this approximation 
is expected to fail for models which include large gradi- 
ents in the thermodynamic variables such as occur at the 
boundaries of shell burning regions. 

The correct identification of the behavior in shell oxy- 
gen burning has wide implications. This stage of massive 
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star evolution is important for a variety of topics of cur- 
rent research interest (e.g., Young, et al. (2005, 2006)). 
In this paper we discuss the quantitative similarities and 
differences between the models of oxygen shell burning 
which we have introduced above and show that the two 
models are in good agreement with each other. In ad- 
dition, the compressible models arc in agreement with 
scaling relations derived from the basic hydrodynamic 
equations as well as analytic solutions to the non-radial 
wave equation for motions in the stable layers. These 
findings lend strong support to the validity of both sim- 
ulations. There arc no signs that either the anelastic or 
the PPM method is breaking down for the conditions 
simulated, even in regions of the flow where the Mach 
number does not exceed M ^ 0.01. 

2. MODEL COMPARISON 

2.1. The Initial Models and Simulation Parameters 

In Table 1 we summarize the initial conditions, com- 
putational domains, zoning, and properties of the devel- 
oped flow for the three models that we will discuss in 
this paper. These include a 2D and a 3D compressible 
model calculated with the PPM method (Meakin & Ar- 
nett 2006a, b) and the non-rotating anelastic model de- 
scribed by Kuhlen, Woosley, & Glatzmaier (2003). The 
initial conditions for the compressible simulations are 
of a 23 M0 star previously evolved with the TYCHO 
stellar evolution code (Young & Arnett 2005) which is 
directly mapped onto the hydrodynamics grid. A 25 
nuclei reaction network is used to track composition 
and energy generation. The computational domain for 
the compressible models are restricted to fractions of a 
sphere and use a spherical coordinate system. The 2D 
model is a 90° wedge embedded in the equatorial plane, 
and the 3D model is a wedge of 30°x30° degrees cen- 
tered on the equator. The radial limits for these models 
enclose both the convectively unstable oxygen burning 
shell, as well as two surrounding stably stratified layers. 
Boundary conditions, which are placed in the stable lay- 
ers, are impermeable and stress free. The net energy 
generation due to nuclear burning and neutrino cooling, 
Lnet = /(ermc + evu)dM « 3.5 X lO"**^ crg/s, is positive 
and goes into PdV work through a background expan- 
sion which develops naturally in the compressible model 
in the same way it does in the initial TYCHO model. 

The anelastic model uses a reference state which is 
a polytropic fit to a 25 M© stellar model which was 
evolved with the KEPLER code (Weaver et al. 1978). 
A multi-fluid model and reaction network are not used 
and nuclear energy generation is instead estimated using 
a power law fit for density and temperature dependence. 
The anelastic model uses spherical harmonics for angu- 
lar coverage, ameliorating the pole singularity problem 
of a spherical coordinate system, and covers a full 47r 
steradians. Chebyshev polynomials are used in the ra- 
dial direction. The radial limits of the anelastic model 
enclose only the convection zone with no regions of sta- 
ble stratification. The boundary conditions, which arc 
within the unstable convectivc layer, arc also imperme- 
able and stress free. The net energy generation is pos- 
itive, L„et ~ 1.5x10^^ erg/s. The anelastic model used 
is unable to model background expansion so the excess 
energy is forced to escape from the outer boundary of 
the calculation. The lower net energy generation of this 



model may be due to it being in a later evolutionary 
stage. 

Both hydrodynamic models use the equation of state 
provided by Timmes & Swesty (2000). The radial limits 
of the oxygen burning convection zone and enclosed mass 
for the two initial models evolved with the KEPLER and 
TYCHO codes are remarkably similar. The compress- 
ible model uses 400 logarithmically spaced radial zones 
(to keep zone aspect ratio dr/rdO ~1) and an angular 
resolution of ~0.3° per zone. The anelastic model uses 
145 zones for a comparable radial extent, and spherical 
harmonics up to order Z = 63 to cover the sphere which 
is roughly equivalent to a Nyquist sampling of '-~^1.5° per 
zone, approximately a factor of five lower angular reso- 
lution than the compressible model. The Rayleigh and 
Reynolds numbers quoted by Kuhlen, Woosley, & Glatz- 
maier (2003) are i?a ~ 5 x lO"^ and Re - 3000. Since the 
compressible model is more strongly driven (larger L„et) 
and has finer zoning (resulting in a lower effective viscos- 
ity), both the effective Rayleigh and Reynolds numbers 
will be higher in the compressible model. 

2.2. Flow Properties: Anelastic Model 

The anelastic simulation has been run for 6500 sec- 
onds. With an average flow velocity of Vc ~ 0.49x10^ 
cm/s, and a radial extent for the convection zone AR f» 
0.39x10^ cm the turnover time tc = 2AR/vc ~ 159 sec- 
onds and the simulation spans approximately 41 convec- 
tive turnovers. The peak velocity is given as Vpeak ~ 
1.8x10^ cm/s which corresponds to a peak Mach num- 
ber M ~ 0.04 for a sound speed Cg « 4.5x10^ cm/s. 
The maximum density fluctuations within the convection 
zone arc found to be of the order p' /{p) ~ 2 x lO^'' which 
is the same order of magnitude as the peak Mach num- 
ber squared, ^ 1.6 x lO^'"* consistent with the scaling 
arguments for the anelastic approximation for thermal 
convection (Cough 1969). 

2.3. Flow Properties: Compressible Models 

The flow in the compressible simulations consists of 
two distinct regimes: the turbulent convection zone, and 
the wave-bearing stably-stratified layers. These can be 
quite readily discerned in the density fluctuation field 
shown in Figure 1. In this section we discuss these two 
regimes in turn for the three dimensional model, and then 
discuss the properties of the two dimensional model. 

The three dimensional compressible model was run for 
^800 seconds of star time. The average flow velocity 
is found to be Vc « 0.8x10^ cm/s. With a convection 
zone width AR « 0.41x10^ cm the turnover time is tc ~ 
103 seconds and the simulation spans approximately 8 
convective turnovers. After an adjustment in the initial 
size of the convection zone due to penetrative convection 
(Meakin & Arnett 2006a,b), the flow achieves a steady 
state within ~200 seconds, or two convective turnovers, 
after which the average flow properties do not change 
appreciably. Figure 2 shows the peak density fluctua- 
tion and the Mach number at each radius for the three 
dimensional model. Within the convection zone, 0.44 < 
r/10^ cm < 0.85, the maximum density fluctuation and 
Mach number are p'/{p) ~ 5 x 10~^ and M ~ 0.09, re- 
spectively. Here also, the fluctuation scale is the same 
order of magnitude as the peak Mach number squared. 
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Fig. 1. — Spatial distribution of density fluctuations are shown for (left) a slice through the 3D compressible model and (right) the 2D 
compressible model. In both models the flow is composed of two distinct regimes, including the convective flow in the center which is 
bounded above and below by stably stratified layers which are host to internal waves. The scale on the 2D model has been truncated to 
same limits as the 3D model for comparison, but extreme values exceed the scale limits by a factor of ~2. The 3D model has been tiled 
twice in angle for clarity. 
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These velocity and fluctuation scales are comparable 
to those of the 3D anelastic model and are listed in Ta- 
ble 1 for both simulations for comparison. The point 
we want to emphasize here is that the character of the 
convective flow is quantitatively in agreement between the 
anelastic and compressible models. We also find, from a 
detailed comparison to be presented separately (Meakin 
& Arnett 2006b), that the 3D compressible model com- 
pares well with the stellar mixing length theory (Kip- 
penhahn & Weigert 1990), including the velocity scale 
{vc ~ lO"^ cm/s) and the superadiabatic stratification 
(AV ~ 5 X 10~ ). The slightly larger velocity scale in 
the compressible model can be attributed to the higher 
Rayleigh number due to the larger luminosity needed to 
be transported by the convective flow. 

2.3.1. Additional Sources of Fluctuations 

In this section we discuss the origins of thermodynamic 
fluctuations which are due to physics not included in 
the anelastic model of Kuhlen, Woosley, & Glatzmaier 
(2003), including stably stratified layers and composition 
effects. We begin with a discussion of the internal wave 
dynamics occurring near the convective boundaries. 

Significant density fluctuations occur at the convective 
boundaries in the compressible model, reaching values 
as large as 11x10^^, over twenty times larger than in 
the body of the convection zone (Figure 2 and Table 1). 
Examining the spatial distribution of the density fluc- 
tuations presented in Figure 1 reveals that the largest 
fluctuations occur at the interface between the convec- 
tion zone and the stably-stratified, wave-bearing layers. 
The morphology of the largest density fluctuations in the 
domain (i.e., those at the convective boundary) are pe- 
riodic in angle and harmonic in time, identifying them 
with internal wave dynamics. In the following discus- 
sion we present analytic estimates for the amplitudes of 
the density fluctuations at the convective boundary and 
show that they are in good agreement with the simula- 
tion data. 

For small amplitude waves the Eulerian density fluctu- 
ations and pressure fluctuations are related by (Unno et 
al. 1989, p.93): 



p' p' 1 N'^ 

—— = — h S^r \- (nonadiabatic terms) (1) 
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with buoyancy frequency N , and Lagrangian displace- 
ment ^r- We defer a discussion of nonadiabatic and com- 
position effects to the end of this section. In the con- 
vection zone, material is nearly neutrally stratified and 
the buoyancy frequency is very close to zero so the sec- 
ond term on the right hand side is not very important. 
At convective boundaries the stellar structure assumes a 
stable stratification with a positive buoyancy frequency, 
and this term can become dominant. This term repre- 
sents the component of the Eulerian density fiuctuation 
due to g-mode oscillations and is the projection of the 
Lagrangian displacements of the wave, onto spherical 
shells. In the presence of steep density gradients, waves 
can lead to large Eulerian fiuctuations even when com- 
pressibility is not important. To be as clear as possible on 
this key point, we give an analogy to well known physics: 
consider waves on a lake. The Lagrangian surface is the 
surface of the water. The Eulerian surface is the aver- 
age level of the water and Eulerian density fluctuations 
occur as the waves (water) and troughs (air) move by 
the observer. The large variation in density is not due to 
compression, but the choice of coordinates (Eulerian in 
this case). 

In order to estimate the amplitude of the density fluc- 
tuations using equation 1 we need to know the size of the 
pressure fluctuation and the maximum radial displace- 
ment amplitude, f""", for wave motions at the bound- 
ary. Both quantities can be estimated by assuming that 
the ram pressure of the convective turbulence is bal- 
anced by wave induced pressure fluctuations at the con- 
vective/stable layer interface: 



'P'w 



(2) 



The validity of this approximation is demonstrated in 
Figure 3, which shows that the RMS horizontal pressure 
fluctuations and the turbulent ram pressure are compa- 
rable in the convection zone and do indeed balance at 
the locations of the convective boundaries. 

The relationship between the pressure fluctuation of 
an internal wave, p^, and the maximum radial displace- 
ment, depends on the wave frequency u and hor- 
izontal wavenumber kh. Perhaps the simplest approxi- 
mation is to assume that internal waves generated at the 
convective boundary are directly related to the convec- 
tion through the convective velocity Vc and eddy scale, 
Ic ^ Hp, by Vc cr/fc/i and kh 27t/Ic in the spirit 
of stellar mixing length theory (Press 1981). Adopting 
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Fig. 2. — (left) Density fluctuations for 3D compressible model: The time average of the maximum density fluctuation is shown as the 
thick line, with the extreme values for the averaging period (two convective turnovers) shown as the shaded region. Within the body of the 
convection zone the average fluctuations are quite low ^^0.25%. The largest fluctuation (the spike at r 0.43 X 10^ cm) is (exceding 
the plot limits), while the secondary maximum (r 0.85 X 10*^ cm) reaches ~3.5%. These fluctuations occur in the radiative regions that 
enclose the convection zone (outside the computational domain of Kuhlen, Woosley, & Glatzmaier (2003)). (right) Macli rmmber for the 
3D compressible model: shown is the instantaneous value of the maximum Mach number at a given radius (thick line), and the rms Mach 
number (thin line). 



these values wo can then use the Unearized momentum 
equation (Unno et al. 1989, p. 96), 



f ma 



h 2 
2^c 



(3) 



and the dispersion relation (for waves in which a <^ N < 
LA, 



kh/kr ~ — 
to estimate the radial displacement: 



^"^^ — ^JJ"*" X kh/kr 



cr2 N 



(4) 



Vc/N. (5) 



The dispersion relation used to connect the horizontal 
and radial displacement amplitudes is valid when the 
wave frequency is much smaller than both the buoyancy 
frequency and the Lamb frequency Li = fc/jCs, with sound 
speed Cs , and is a reasonable approximation for the wave 
properties adopted above (where a/LiK, Mc). 

Here, our main result is the last expression in equation 
5 for the radial displacement amplitude of the interfacial 
wave which is in pressure balance with the ram pressure 
of the convection. This expression is equivalent to the 
statement that the kinetic energy of the turbulent motion 
exciting the wave is balanced by the potential energy of 
the wave, which follows naturally from the basic energetic 
properties of waves in fluids (Lighthill 1978). 

Finally, we use the displacement given by equation 5 
and the pressure fluctuation in equation 2 with equation 
1 to arrive at our estimate for the interfacial density fluc- 
tuation amplitude, 



VcN 



(6) 



in terms of the Mach number Mc, gravity g, and buoy- 
ancy frequency N. Adopting flow parameters from the 

simulation (?;,, ^ 10^, g ~ 10^ in cgs units) we find 
p'/{p) ~ (10^3 + VV X 10^2) ^ith N in rad/s. The va- 
lidity of this expression is apparent when comparing the 



buoyancy frequency in Figure 3 with the density fluctu- 
ations at the convective boundaries in Figure 2. 

It can also be seen that the density fluctuations 
throughout the stable regions, not just at the convec- 
tive boundaries, are in rough agreement with the scal- 
ing given by equation 6, though the amplitudes in the 
stable layers drop off with distance from the convective 
boundary. This is due to a spectrum of internal wave 
modes excited at the convective boundaries, rather than 
the single mode assumed in the above analysis. Each 
modal component contributes to the pressure balance at 
the convective boundary and is composed of wave packets 
that travel back and forth within the resonating cavity 
of the stable layer (Unno et al. 1989) causing fluctua- 
tions throughout the region. A more detailed analysis is 
possible in which the entire spectrum of internal waves 
is estimated by matching the wave motions to those of 
the spectrum of turbulent convection (e.g. Carruthers & 
Hunt 1986). Our single mode approach, however, works 
very well in describing the amplitudes of the fluctuations 
at the convective boundaries and provides a reasonable 
upper limit to the fluctuations throughout the entire sta- 
ble layer. 

We conclude this section with a discussion of the role 
that entropy fluctuations play in setting the scale of den- 
sity fluctuations and hence buoyancy of material in the 
convection zone. The non-adiabatic term in equation 1 
takes the form, 



Pi 

iP) 



Vt 
Cp 



6S 



(7) 



with thermodynamic derivative Vt = —{din p/dlnT)p, 
specific heat at constant pressure Cp, and Lagrangian en- 
tropy fluctuation 5S. In the present model the largest 

non-adiabaticity is due to net effect of nuclear burning 
and neutrino cooling. The entropy fluctuation can then 
be written SS„uc = ^Qnei/T ~ enei<^t/T where the time 
material dwells in the burning region is St « Ar/vc ~ 
Ic/Vc ~ 10s. For the nuclear energy release in the cur- 
rent model, with a peak burning rate of e„et ~ 10^** 
erg/g/s, the maximum entropy fluctuation will be of or- 
der SSl^^^c ~ 0-01 ill units of NAks- The correspond- 
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Fk;. 3. (left) I'rcssure fluctuations in 3D compressible model: The time averaged horizontal RMS pressure fluctuations are shown as 
the thick line, with the envelope of extreme values over two convective turnovers indicated by the shaded region. The radial ram pressure 
of the turbulent convection, pv'^, is shown as the thin line. The curves cross at the convective boundaries where the turbulent pressure is 
balance by the pressure fluctuation induced by internal waves in the adjoining stably stratifled layers, (right) The buoyancy frequency is 
shown in units of rad/s. Also shown as the dashed line is the buoyancy frequency normalized by the gravity which, through Equation 6, 
sets the scale of the density fluctuations at the convective boundaries (compare with Figure 2a). 



ing maximum density fluctuation will be approximately 
p'/{p)\nuc ^ 2.5 X lO^'', which is of order the amplitude 
of the fluctuations due to the weak compressibility ef- 
fects in the convective flow. These values compare well 
to the perturbations estimated using mixing length the- 
ory, consistent with the convective flow being driven by 
the nuclear burning luminosity in the shell. 

Entropy fluctuations also occur within the convection 
zone due to composition inhomogeneities. The radial en- 
tropy profile and the RMS entropy fiuctuations are pre- 
sented in Figure 4. The fluctuations are due to: (1) inter- 
facial wave motions which cause Eulerian fluctuations in 
the same manner as for the density fluctuations discussed 
above; and (2) the entrainment of high and low entropy 
material at the convective boundaries which is mixed into 
the nearly adiabatic convection zone. The wave induced 
flucttiations appear as spikes near the convective bound- 
aries and are present in both the curve of minimum and 
maximum fluctuation. The regions affected by composi- 
tional inhomogeneities are labeled in Figure 4, with low 
entropy material entrained from below and high entropy 
material entrained from above. The entropy fluctuations 
associated with this material is another source of density 
fluctuations and explains the larger values that occur just 
within the boundaries of the convection zone in Figure 2. 
The entropy fluctuations associated with the entrained 
material are much larger than due to nuclear burning 
(the entropy perturbation in the region of greatest nu- 
clear energy deposition, r ~ 0.45 x 10^ cm, is primarily 
negative). The entrainment of material from stable layers 
by a turbulent convective flow is an essential addition to 
stellar evolution modeling with signiflcant consequences 
for the evolution of burning shells in presupernova mod- 
els. 

2.3.2. Dimensionality 

It has long been known that 2D simulations were in- 
formative only to the extent that care is used in their 
interpretation. In 2D the vorticity is restricted to the 
direction normal to the computational domain, while in 
3D instabilities cause its orientation to wander. Thus 
2D is useful in situations in which there are physical rea- 
sons to enforce the symmetry (e.g., terrestrial cyclonic 
storms), but has nevertheless been used widely in more 



general applications because of computer resource limi- 
tations. The increasing availability of computing clusters 
and software parallelization tools is now making 3D hy- 
drodynamic simulation more common, and we are start- 
ing to assess the adequacy and limitations of earlier 2D 
work. 

We have calculated a 2D compressible model for 2400 
seconds of star time. We flnd an average flow velocity 
of Vc ~ 2.0x10'' cm/s and a convective turnover time of 
tc ~ 40 seconds, so our simulation spans approximately 
60 turnover times. The peak velocity during the course 
of the simulation is '-^^5.5x10'' cm/s, corresponding to a 
peak Mach number of M ^ 0.163. The density fluc- 
tuations within the convection zone reach a maximum 
value of p^/(p) ~ 6xl0~^. At the convective boundaries 
the density fluctuations attain a peak value of pj,/ (p) ~ 
12x10-2. 

We find two significant differences between the 2D and 
3D models. First, we find a significantly decreased tur- 
bulent mixing rate in the 2D simulation. Material en- 
trained into the convection zone at the boundaries is 
pulled into the large cyclonic flow patterns in the 2D sim- 
ulation where large composition inhomogeneities persists 
for several convective turnovers. In contrast, material en- 
trained into the convection zone in the 3D models is ho- 
mogenized within a single convective turnover time. This 
effect is illustrated in Figure 5, which shows the spatial 
distribution of oxygen abundance, as well as RMS fluc- 
tuations for both the 2D and 3D simulations. The 2D 
simulation retains high level fluctuations throughout the 
convective zone, while inhomogeneities in the 3D model 
are mixed to low levels by the time material completes 
a single crossing. Comparing Figures 5 and 1, which 
are snapshots at the same time, show the high entropy 
oxygen entrained at the top boundary corresponds to a 
negative density perturbation. 

The second major discrepancy between the 2D and 3D 
models is the convective velocity scale. We find that 
both the mass averaged convective velocity and the peak 
velocity fluctuation are ~ 2 time larger in the 2D model. 
This velocity scale difference may be connected to the 
lower ttirbulent mixing efficiency in the 2D flow. We 
find that the net enthalpy flux, which consists of upward 
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Fig. 4. — (top) Entropy profile of for 3D compressible model, (bottom) Entropy fluctuations for 3D compressible model: The two solid 
lines indicate the maximum and minimum fluctuation for a given radius over the course of two convective turnovers. The annotations 
indicate fluctuations due to low entropy material entrained at the lower boundary and high entropy material entrained at the upper 
boundary. 



and downward directed components, F„et = Fup — Fdown, 
is the same between the 2D and 3D models. In the 2D 
model, however, the relative value of the individual flux 

components relative to the net flux, e.g., Fy_p/Fnet, are 
much larger than in the 3D model. Therefore, the 2D 
model requires a larger velocity scale to move the same 
net flux due to the inefficiency of depositing the advected 
enthalpy across the convection zone. 

We find that 2D and 3D models compare well in the 
wave region, but differ in the convection zone. The 3D 
convection is more similar to that of the anelastic model 
of Kuhlen, Woosley, & Glatzmaier (2003) and the val- 
ues predicted by mixing length theory. While waves be- 
have similarly in 2D and 3D, turbulent convection doc^s 
not, particularly with regards to turbulent mixing effi- 
ciency. Although the spatial resolution of the 2D and 
3D models are the same, the number of degrees of free- 
dom in the angular direction is much larger in the 3D 
model, N3D/N2D = (100)7320 -31. If the muiiber of 
degrees of freedom were the only important parameter, 
one might wonder if 2D would provide a more efficient 
surrogate to 3D. It turns out that 2D is actually more ex- 
pensive than 3D because the same degrees of freedom in 
2D requires a higher spatial resolution and hence a more 
severe time step constraint, and computational cost is 

^cost N space ^ ^time' 

3. CONCLUSIONS 

A comparison between the flow properties in fully com- 
pressible and anelastic simulations of stellar oxygen shell 
burning indicate that the two methods produce quantita- 
tively similar results. Both methods produce convective 



flows in 3D models which are compatible with the results 
expected for the mixing length theory of convection for 
this phase. The compressible models have been extended 
to include additional physics not included in the anelas- 
tic model, namely stably stratified boundary layers and a 
multi-fluid flow {Ngpedes = 25). The interaction between 
the convection and the stable layers excite internal waves 
which produce larger thermodynamic fluctuations (up to 
11% in 3D). Composition inhomogeneities due to ongo- 
ing entrainment events at the convective boundaries also 
cause density fluctuations on the several percent level, 
though material is homogenized rapidly in the 3D model 
through turbulent mixing. 

The relatively large fluctuation which arise at the con- 
vective boundaries ~ 11% may stress the reliability of the 
anelastic approximation if this region is to be included 
in future simulations of oxygen burning or later epochs, 
where entropy and density gradients are large. A variety 
of convection studies have shown that boundary condi- 
tion type (e.g., hard wall compared to stable layer) alters 
the overall flow pattern within a convection zone (Hos- 
sain & MuUan 1993; Rogers & Glatzmaier 2005b) and 
therefore the astrophysically correct conditions should 
be used. Low Mach number solvers (e.g. Lin et al. 2006) 
may be the most efficient tools for extending studies of 
oxygen and silicon shell burning to full spherical domains 
in 3D while retaining the crucial density gradients at the 
convective boundaries where convective penetration and 
entrainment operate, and asymmetric fluctuations arise 
which may have important implications for the evolution 
of pre-supernova models. Earlier stages such as carbon 
and neon burning have both milder flows and shallower 
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Fig. 5. — Spatial distribution of the oxygen mass fraction is shown for (left panels) the 3D compressible model and (right panels) the 
2D compressible model. The spatial distribution is shown in the top row, and the time averaged RMS horizontal fluctuations are shown in 
the bottom row wih the shaded region indicating extreme values of the fluctuations over two convective turnovers. 



density gradients and should be better suited for anelas- 
tic methods, even at convective boundaries. Background 
expansion and multi-fluid effects should be included how- 
ever. The large time step advantage of the anelastic and 
low Mach number simulations allows for much larger do- 
mains or better resolution. 

Although the efhciency of fully compressible hydrody- 
namics may be low for the Mach numbers modeled, there 
are no signs that the solver used is breaking down in the 
oxygen shell burning simulations presented here. This 
conclusion is supported on several grounds, including: 
(1) The compressible model is in good agreement quan- 
titatively with the anelastic methods for the convection 
zone region, including the velocity scales, and thermo- 
dynamic fluctuation amplitudes, a region in which the 
anelastic method is expected to perform well. (2) The 
compressible simulation of the convection zone is also in 
good agreement with the results of the one-dimensional 
TYCHO model, including the velocity scale and back- 
ground stratification estimated using mixing length the- 
ory. (3) The dynamics in the stably stratified layers 
in the simulation agree well with the analytic solutions 
to the non-radial wave equation, including the decom- 
position of the flow into speciflc, unambiguous modes 
(Meakin & Arnett 2006a). (4) The fluctuation ampli- 
tudes at the convective boundaries which are due to wave 
motions are found to be in good agreement with analytic 



estimates for their scale. 

Contrary to the assertion made by Almgren et al. 
(2006) that compressible codes should fail for M < 10~^, 
we find a robust solution that agrees with an anelastic 
method for the same region simulated. Additionally, re- 
cent compressible simulations of He shell fiasli convec- 
tion by Herwig et al. (2006) using the finite-volume Go- 
dunov code RAGE (Baltrusaitis et al. 1996) find a flow 
with M ~ 10"'^, with apparently robust results, includ- 
ing well behaved g-modes. Almgren et al. (2006) present 
an example simulation illustrating the failure of PPM to 
track temperature for a simple flow with Mach number 
Af ^ 0.05. This calculation, however, uses a compress- 
ible PPM code (FLASH) with two major differences from 
ours (PROMPI): they used the hydrodynamic procedure 
described in Zingale, et al. (2002) to remove the hydro- 
static pressure from the Riemann solver, and their stellar 
model was much more degenerate than ours (the equa- 
tion of state tends to become independent of tempera- 
ture under their conditions, and care must be taken with 
cancellation of terms). 



This work was supported in part by the ASCII FLASH 
center at the University of Chicago. One of us (DA) 
wishes to thank the Aspen Center for Physics for their 
hospitality. 
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TABLE 1 

Comparison of Oxygen Burning Models 



Variable 


Units 


2D-PPM 


3D-PPM 


SD-Anelastic'^ 




(M© ) 


23 


23 


25 




(M© ) 


1.0(1.5), 2.7(2.4) 


1.0(1.5), 2.7(2.4) 


1.2, 2.3 




(10'' cm) 


0.3(0.44), 1.0(0.85) 


0.3(0.44), 1.0(0.85) 


0.45, 0.84 


Lnet^ 


(crg/s) 


3.2x10'*'^ 


3.5xl04« 


1.5x10*6 


AO,A<l> 


(dcg.) 


90 


30, 30 


360, 180 


Zones/Modes 


(rirXn^xng) 


400x320x1 


400x100x100 


145x63(0 x31(m) 




(lO'^ cm/s) 


7.2 


3.8 


1.8 


{Vconv ) 


(10'^ cm/s) 


2.0 


0.8 


0.49 






0.163 


0.09 


~0.04 






0.03 


0.01 


~0.01 


tconv 


(s) 


40 


103 


159 


tmax 


(s) 


2400 


800 


6500 


max{ p'Jip) }d 
max{ Z/{T} } 


10-2 


1.0(6.0)"^ 


0.5 


0.2 


10-2 


0.25 


0.05 


0.06 


P'b/{p) } 
max{ Tj'/(T) } 


10-2 


12.0 


11.0 




10-2 


2.7 


1.0 





^The values quoted for the 3D anelastic model are from Kuhlen, Woosley, & Glatzmaier (2003). 

^Zero age main sequence mass. 

^Values in parentheses indicate the extents of the convection zone for the compressible models and 
the other values indicate the extents of the entire computation domain including the stable layers. 

^The net luminosity for the compressible models, I-inet = fi^nuc + €t^p)dM, is estimated at t~400 s 
and is slowly decreasing with time due to an overall background expansion occuring within the burning 
region. 

'^The c and b subscripts indicate thermodynamic fluctuation amplitudes that are estimated in the 
convection zone and in the region of the convective boundary, respectively. 

■^Thc value in parentheses is the maximum density fluctuation estimated over two convective turnovers 
while the other value represents the time averaged maximum fluctuation. The significantly larger 
fiuctuations seen in the 2D model compared with those in the 3D model occur near the center of large 
vortices which persist for several convective turnover times in the 2D model but are absent in the 3D 
convective flow. 



